Josephson current through a molecular transistor in a dissipative environment 
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We study the Josephson coupling between two superconductors through a single correlated molec- 
ular level, including Coulomb interaction on the level and coupling to a bosonic environment. All 
calculations are done to the lowest, i.e., the fourth, order in the tunneling coupling and we find 
a suppression of the supercurrent due to the combined effect of the Coulomb interaction and the 
coupling to environmental degrees of freedom. Both analytic and numerical results are presented. 
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I. INTRODUCTION 

Mesoscopic systems connected to superconducting 
leads have been investigated for a number of years. If the 
mesoscopic system itself is superconducting the transport 
is influenced by the so-called parity effect in the Coulomb 
blockade regime*^ For normal metal grains, the effect of 
the superconductivity is merely to modify the tunnel- 
ing density of states due to the density of the supercon- 
ducting leads, thus introducing an additional gap in the 
current- voltage characteristics. 

Naturally, the Josephson current is also affected by the 
Coulomb correlations, because the transfer of a Cooper 
pair charges the mesoscopic system by two electron 
charges. This process is similar to a cotunneling event, 
and in the limit of large Coulomb repulsion, the path that 
involves double occupancy of the central region is not al- 
lowed, which results in a suppression of the Josephson 
coupling. This has been studied extensively in a num- 
ber of theoretical works, starting with the work of Shiba 
and Soda 3 " and Glazman and Matveevpi who calculated 
the Josephson current in the limit of infinite Coulomb 
repulsion and in perturbation theory, the leading order 
being the fourth order in the tunneling amplitudes. Inter- 
estingly, the Josephson current was shown to change its 
sign when the level which supports the tunneling current 
becomes occupied. This so-called 7r junction behavior, 
which is consequently also relevant in the Kondo regime, 
has been studied in a number of papers^ii For larger 
dots the interplay of the above mentioned parity effect 
and the Josephson effect has also been addressed^ 

Experimentally, there has been a number of stud- 
ies of metallic wires connected to superconductorsji 3 ^ 
but only a few studies in the Coulomb blockade regime. 
Buitelaar et ulJ±H± have observed multiple Andreev re- 
flections and Kondo physics in carbon nanotube quantum 
dots. The supercurrent was not directly observed. 

In this paper, we study the Josephson coupling through 
a single level system coupled to vibrational or dissipative 
environments. This is relevant for molecular transistor 
systems where strong influence of the coupling to various 



vibrational modes have been observed m2r— In the case of 
normal metal leads this has lead to a number of related 
theoretical worksi^Sr— However, so far the combination 
of the vibrational coupling and correlation transport has 
not been studied in the context of supercurrent. Be- 
ing a groundstate property the dissipative environment 
is expected to have a more dramatic effect on the super- 
current as compared to usual electron tunneling. This 
is indeed what we find in the case of strong coupling, 
where the Franck-Condon factors strongly suppress the 
supercurrent. 

The paper is organized as follows. In the next sec- 
tion our model Hamiltonian is presented and a convenient 
unitary transformation is performed. In Sec. Ill, the ba- 
sic formula for the Josephson current is derived to the 
fourth order in the tunneling amplitudes, and the differ- 
ent ingredients in this formula are discussed. Section IV 
gives the supercurrent without coupling to the bosons, 
both with and without Coulomb interaction, while the 
Sec. V studies the full case and discusses different limit- 
ing cases. Finally, in Sec. VI our conclusions are stated. 
The technical details of analytic calculations are put into 
two Appendixes. 



II. MODEL HAMILTONIAN 

The model we study in this paper consists of two su- 
perconducting leads and a central region described by 
a single level with Coulomb interaction and its charge 
occupation coupled to one or many harmonic oscillator 
modes, see Fig. ^ 

The model Hamiltonian is given by 

H = Hq + H T , (1) 
where the unperturbed part of the Hamiltonian is 

H = H L + H R + H M , (2) 
with Hl r being the BCS Hamiltonians for the left and 
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where 



FIG. 1: Schematic picture of the physical setup described by 
the model. The central region (molecule/quantum dot) is cou- 
pled to two superconducting leads and can be independently 
gated by a gate. Electrons entering the central region expe- 
rience mutual Coulomb interaction and interact with oscilla- 
tory modes — either one distinguished mode such as center- 
of-mass vibration of the molecule as a whole or many modes 
due to intramolecular vibrations or substrate phonons acting 
as a heat bath. 



right leads, respectively, 

H a = E e ^c\ aa c kaa - E A « c Lt c -fc<4 + H - c - ' 

(3) 

with a = L, i?, and the complex gap functions defined as 
A Q = lAoJe^". Furthermore, Hm describes the central 
region 



Hm = H 



MO 



Hvih ~\~ H e — vib; 



(4) 



where Hmo is the Hamiltonian for the electronic degrees 
of freedom 



X 



(8) 



In this model we neglect any modification of the tunnel- 
ing due to the vibrations (this is very often the relevant 
situation and the model could easily be generalized to in- 
corporate such a dependence) , and the tunneling Hamil- 
tonian is therefore given by 



a=L,R 



(9) 



where Hi 

1 a 



H, 



To) 



and 



^T Q -E ff T QCT > ^t^-E^L^- (io) 

<r=T4 k 

In the next section we will calculate the Josephson cou- 
pling using perturbation theory in the tunneling. For this 
purpose, it is convenient first to use a polaron represen- 
tation, which transforms the coupling term i? e _ V ib to a 
displacement operator in the tunneling termi2&2& 

The unitary transformation 



H = SHS\ S = e- tA{n t +n i\ 



where 



A — ^ ^ IjyPty, i v — 



(11) 



(12) 



removes the coupling term _ff e _ v ;b from the Hamiltonian 
at the expense that the tunneling term acquires an os- 
cillator displacement operator, so that Eq. lfTU|l becomes 



Hmo = £ E d ° dcr 



(5) 



H T a a = = E* fef * C W d ° 



(13) 



with d(j , d\ being the operators for the local level, and n-\ 
and nj_ the corresponding spin-dependent occupations. 
We assume that the single-particle energy level £ can be 
experimentally tuned by the gate, i.e., £ = £(V g ), and 
present the results for the critical Josephson current as 
functions of £. We have here defined the origin (zero) of £ 
to be at the Fermi level of the superconducting electrodes, 
i/vib is the Hamiltonian for the vibrational degrees of 
freedom 



H, 



vib 



E 



p 



2m, 



(6) 

and, finally, i? e _ v ib describes the coupling between elec- 
tronic and vibrational degrees of freedom. This coupling 
is solely through the charge on the central system, i.e., 



(n T +n l )X, 



(7) 



Furthermore, the transformation renormalizes the on-site 
energy and the Coulomb interaction according to 



!=£-^E A ^' i> = c/-E A ^- 



(14) 



In the following we will skip the tildes and use just £, U 
again but we mean the renormalized quantities. Using 
this transformed Hamiltonian, we calculate the Joseph- 
son current to the lowest order in the tunneling Hamil- 
tonian in the following sections. 



III. JOSEPHSON CURRENT 

The current operator for the current through contact 
a = L, R is N a (we use K = e = 1 throughout the whole 
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paper), where N a is the operator of the number of elec- 
trons in lead alpha. After the unitary polaron transfor- 
mation introduced previously, we hence obtain the cur- 
rent, I a , as 



In, = i 



H,N Q 



= i(H+ a -H^ a )= 21m (H Tc 



(15) 

Performing the standard thermodynamic perturbation 
expansionist in the tunneling we obtain for the Joseph- 
son current in the lowest non-vanishing order, which is 
the fourth order, in 



1 f 13 f 
7 Q = -2Im — / dn / dr 2 I dr 3 
•J' Jo Jo Jo 

T T (H T ( Tl )H T (T 2 )H T {Tz)H Ta 



P 



(16) 



The Josephson current must involve two an d two 
H^p , which can be chosen in three ways, and hence 



rP rP rP 

I a = -Im / dr\ I dr 2 \ dr 3 
Jo Jo Jo 



X 



Tr m & (rim & (.r 2 )H^(r 3 )H Ta 



(17) 



where we also used that in order to have Cooper pair 
tunneling, the must belong to the junction opposite 
to where the Josephson current is "measured" via Hj^ a , 
i.e., a means the lead opposite to a. Because of spin 
symmetry we can choose the spin of the last as, say, 
spin up, which then means that the other carries spin 
down. In the same way, the spin of the two can be 
chosen in two ways. All in all, we thus obtain 

rP rP rP 

I a = -41m / dri j dr 2 / dr 3 
Jo Jo Jo 

X 



rP rP rP 

= -41m tl & t*_ p& t-k a tka / dn / dr 2 / dr 3 

<. „ Jo Jo Jo 



k p 



x Fk a {T3)F* & {Ti - t 2 )B{t x .,t 2 ,t 3 )V{tx,t 2i t 3 ), 

(18) 

where T arc the anomalous Green's functions of the leads 



= - (Tr (cL fcai (r)4 aT (0))) o (19) 

and where we define the following two functions pertain- 
ing to the central region: 

8(7-1,7*, r 3 ) = (T T (4(r 1 )d t T (r 2 )rf i (r 3 )d T (0))) o (20) 



and 



^(ri,T 2 , r 3 ) - (T T {e- iA ^e- iA ^e iA ^e iA )) 



The first function B describes the propagation of a 
Cooper pair through the central region, while the other 
function T> accounts for the corresponding shifts of the 
oscillator degrees of freedom when the charge on the cen- 
tral region is changed. 

The anomalous Green's functions T are easily calcu- 
lated in the standard way^SiSl and we have 



Fka(T) = - " fa{Eka,T), 
Ltjku 



(22) 



where 



f a (E ka ,r) = e - E *<*M - 2 coshtEteT) n F {E ka ) (23) 



and as usual E ka — \/ z\ a + |A Q | 2 . Throughout, we will 
assume low temperatures such that | A^.i? |/? 1, and we 
can thus approximate 



fa(Eka,T) S3 e 



-E ka \r\ 



-E ka (P-\r\) 



(24) 



In expression (|18fl for the current there are two sums 
over states in the superconductors, which define the tun- 
neling density of states (TDOS) as 

r a (e) = 2n^2t ka t^ ka 5(£~E ka ) = 2tt \t ka \ 2 S(e-e ka ). 

k k 

(25) 

For a small central region, where the coupling is point- 
like, we can approximate t k by a constant, which gives a 
weak energy dependence of Y. 



A. Critical current 

The Josephson current is a function of the phase- 
difference between the two superconductors and using 
Eqs. (EHJ, <l22ll . and (55j), we have 

I a = I c sin 6, (26) 

where the phase difference 8 = a — (f> a . Finally, the 
critical current I c is given by 



TT 2 , 



00 poo 



1 f 13 r 13 [ p 

dri / dr 2 / dr 3 / de j de' 



—00 j — 00 



(21) 



x T L (e)T R (s') ] -^^f L (E,T 3 )f R (E',n - r 2 ) 
x B(n ,t 2 ,t 3 ) V{ti , r 2 , 73). 

(27) 

This expression forms the basis for the further calcula- 
tions in this paper. In fact, from now on we will assume 
that the two tunneling densities of states are energy in- 
dependent. 

The value of I c may come out negative and we will see 
in the following that it really does so due to the Coulomb 
interaction. The case with I c < is called a it junction. 
This terminology originates from an equivalent descrip- 
tion of the Josephson junction in terms of total energy 
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(or free energy at nonzero temperature) of the junction 
as a function of the phase difference E(9). Since I a — 
2dE(9)/d6 the total energy reads E(9) = -(/ c /2)cos6> 
and reaches minimum at 9 = for I c > or 9 = it fol- 
ic < 0, respectively (i.e., the ground state of the junction 
corresponds to equal/opposite phases in the two leads). 
The 7r junction behavior has been noted in a number of 
papers^ii The origin of this sign change is the blocking 
of channels for the Cooper pair exchange when U is large. 
See Ref. for a detailed account. 

Ideally the Josephson current can be measured in a 
current-bias setup where there is no voltage drop across 
the junction until the critical current is reached. How- 
ever, in practice current-bias is difficult to achieve for 
large resistance junctions such as these single-electron 
devices and instead a voltage is swept across the junc- 
tion and the critical current must be determined as half 
of the area of the (ideally ^-function-like) peak in the 
dl / dVcvxve around V — 0, which then, of course, only 
yields the absolute value of I c . 



B. Function B 

Since the central region has interactions, we cannot 
use Wick's theorem, and we must evaluate the function 
B using the many-body states, of which there are four: 
|0),|T),U), and |U). In Eq. CO only ||) and |U) con- 
tribute to the trace, i.e., B = B\ + B2, where 



TV (4(^)4(7-2)^(7-3)^(0)) I ?) o , (28a 



B-2 = P U (t! \Tr (4(71)4(72)^(73)^(0) 



with 



1 + 2e-« + ' 



P 



T[ 



1 + 2e-K + e -pE 2 ' 



(28b) 

(29a) 
(29b) 



and E2 = 2£ + U. For B\ only three orderings of the 
operators give a nonzero result, and we find 

Bi = P T {e^- T3+ ^9(T 2 - t 3 )9(t 3 - 7i) 

+ e i32(r 2 -T 3 ) e e(r 1 -r 2 +r 3 ) e(T3 _ _ n) 

+ e B2(r 1 -T3) e C(r 2 -r 1 +r 3 ) 6) ( T3 _ ^Qfa _ T2 )| 



(30) 



Likewise, for B2 we find 

B 2 = -P U {e^- T °- T ^e E ^6(n - t 2 )9(t2 - r 3 ) 
+ e^~ T2 - T3 *>e E ^9(T 2 - Ti)6(n - 73) 



(31) 



C. Function V 

We also calculate the function T> in Eq. (|21l) involving 
the bosonic degrees of freedom. Since the unperturbed 
Hamiltonian is quadratic in boson operators, this is eval- 
uated as^i 



(32) 



V = exp [h{ Tl - 73) + h(n) + h(r 2 - t 3 ) + h(r 2 ) 

- h(n - t 2 ) - /i(t 3 )Y 

where the function h is defined as 

h(r) = (T T (A(r)A(0))) -(A 2 ) . 



(33) 



When inserting the formula l|12|) for the operator A one 
obtains 



h(r) 



dio 



(n B H(e u|r| -l) 



+(l+n s (o;))( e -H_i)) 
with the spectral function of the bath J(u>) 



(34) 



2m u io v 

and the Bose function «s(w) 
n B {u>) 



1 



exp(/3w) — 1 ' 



(35) 



(36) 



In this paper, we concentrate on two cases — either 
a single important vibrational mode or a bath of har- 
monic modes. For the latter case, we consider mainly 
the situation where the dispersion relation corresponds to 
the so-called Ohmic case, equivalent to a frequency inde- 
pendent damping coefficient or, equally, spectral function 
J (to) linear in frequency. 



1. Single oscillator case 

The case of a single oscillator with frequency luq , mass 
m and coupling constant A corresponds to J osc (o;) = 
gu>QS(uj — wo), where the dimensionless g = A 2 /2mo>g as 
in Ref. (23J, and we have 

/W(r) =g[n (e^ - 1) + (1 + n )( e -^ - 1)] , (37) 

with uq = ng(wo). 

2. Ohmic environment case 

The other generic case that we consider is that of 
Ohmic heat bath in which case we use the spectral func- 
tion 



JohmM = gujexp(-uj/uj c ) 



(38) 
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parametrized by the dimensionless interaction constant 
g and the upper cutoff frequency u> c . In this important 
case it is possible to calculate at T = the correlation 
function h(r) from Eq. (|34|) analytically yielding 



and 



h hm(T) = -g\n(l + w c \t\ 



(39) 



IV. 



JOSEPHSON CURRENT WITHOUT 
DISSIPATION 



As a reference for the discussion on the influence of 
coupling to a dissipative environment, we discuss the case 
without coupling to the bosonic environment. This we do 
in three steps: first we set U = 0, then we look at the 
infinite U case and finally we give the expression for the 
general case. 



A. No Coulomb interaction (7 = 

This result is derived in Appendix and the critical 
current is found to be 



2(1 



T L T R / tanh(/%/2) 
- (C/A) 2 ) V £ 



tanh(/3A/2) \ 
A /' 
(40) 

We note that this expression diverges when £ = and 
T = 0, which however is regularized by higher order 
terms in T's, see discussion in Appendix [S] Further- 
more, we see that the critical current is always positive, 
which means that the negative critical currents, i.e., the 
7r junction behavior, found below are a result of the cor- 
relations. In Fig.[5]this noninteracting result is compared 
with two interacting cases U = A and U — > oo. 



B. Strong Coulomb interaction U — ► oo 

Let U — > oo such that the doubly occupied state is 
taken out. Hence 

Bi = P T e« (Tl - T3+r2) #(T 2 - t 3 )9(t 3 - n ) , (41) 

and when inserting this into formula for the critical cur- 
rent Ij27(l . performing the imaginary time integrations 
using the approximation (|24|l and ftA ^> 1 (so that 
exp(-ftE ka ) « 0), we find 4 ^ 



(42) 



where 
AT(0 = 



|A L |d£ 



\A R \dE> 



Al \ ^/E 2 -\A L \ 2 J\ Ar \ VE' 2 - \A R \ 2 



:C{E,E') 

(43) 



c(£,£')= T 



2e~« / e« 

1 



(44) 



At T = 0, this reduces to 



f -2 + + £)(£' + £)] fore>0, 

(45) 

The integral in Eq. (|43|) can be performed analytically at 
T = 0, and |A L | = |A fl | = A yielding 

f 2n(£/A)/A for £ > 0, 

JV(0 = S ( 46 ) 
\ -n(-£/A)/A for £ < 0. 

The dimensionless function n(x) defined as (x > —1) 



n{x) - 



du 



dv 



1 



/i Vu 2 - 1 J i Vv 2 - 1 (u + v)(u +x)(y + xy 

(47) 

can be expressed by (see Appendix lB|l 

^-(1 — x) — arccos 2 x , . 

n(x)= 41 x{ l_ x2) . with -Kx, (48) 

where the analytic continuations of arccos x = i ln(a; + 
V 'x 2 — 1) for x > 1 is understood. The function n(x) 
is always positive, it diverges at x —> — 1 + , and then 
smoothly decays for increasing x with the asymptote 
n(x) ~ 7r 2 /4x 2 for x — > oo. The expression (|46|) is 
compared with the noninteracting and finite U results 
in Fig. [21 We see that the magnitude of the critical cur- 
rent is highly suppressed by the very strong Coulomb 
interaction. 

For finite but small temperatures ft A > 10 so that our 
approximations are still valid, we can write for the magni- 
tude of the Josephson current (for |£| < A only, otherwise 
the zero temperature expression should be used) 



2r L r fl , n(£/A) - e"«n(-£/A) 



7T 2 A 



(l + 2e-«) 



(49) 



The finite temperature behavior in the limit U — > oo is 
illustrated in more detail in Fig. [3] for three values of 
temperature. In the lowest temperature curve ft A = 10 
we also compare the analytic expression l|49[) with a direct 
numerical evaluation of the triple imaginary-time integral 
(Eq. H53fl with V = 1) routinely used for the dissipative 
cases. We see an excellent agreement between the two 
methods. 

In an experiment with a single Josephson junction the 
absolute values of the presented curves would be mea- 
sured. This would give curves with a dip down to zero 
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at £ — (for a finite temperature) and with asymmet- 
ric shoulders around the dip with the ratio between the 
shoulder heights being 2. Even though the dip may be 
smeared in the experiment for low enough temperatures 
the asymmetric shoulders should persist thus revealing 
the crossover to the tt junction regime. 



C. Finite Coulomb interaction U 

For a finite value of U we have to consider all terms of 
B = B\ + £>2 in an analogous way as previously the first 
one in the U — > oo case. Neglecting again terms e~^ Ek - p , 
we recover the formula (|42|l for the current, but with the 
function C(E, E') in Eq. lO replaced by 



C{E,E') = 



Z(E + E') 



L(E,E r ), 



(50) 



where Z=(l + 2e~K + e ~mt+u)) and 
L(E,E') = - 
+ 



-0i 



(51) 



(E + 0(E' + {E - t){E> - 

2e-^ +u ^E + E l ) 
(2t + U)(t + U-E)(£ + U-E>) 

2{E + E') 
(2Z + U)(Z + E)(Z + E>) 

2e-K 
[E'-Z)[E' + Z + U) 

(E-Z)(E + Z + U) 



{E + £ + U){E' + £ + U) 

W+U -E){£ + U -E'Y 



The resulting integral does not seem to be analytically 
calculable in the whole £-range, not even at T = and 
= | A,r | = A. Yet, in that limit one can achieve 
significant simplifications at least for some £'s which even 
allow us to evaluate the n(x) function defined by Eq. H47|) 
yielding Eq. 1)48 [I. The details of those calculations can 
be found in Appendix El We have, however, calculated 
the critical current numerically, and an example is shown 
in Fig. H for U = A. 



JOSEPHSON CURRENT WITH 
DISSIPATION 



Next, we study how the critical Josephson current 
changes, when coupling to vibrational modes is included. 
In order to do that, we will perform a numerical inte- 
gration of the three imaginary time integrals. For this 




FIG. 2: (Color online) Josephson current dependence on the 
Coulomb interaction strength at zero temperature. Shown 
are the critical Josephson current through an Anderson level 
with no Coulomb interaction U = (full line), with moder- 
ate interaction U = A (dash-dotted line), and with infinite 
repulsion U — > oo (dashes). The interacting cases exhibit the 
phenomenon of the tt junction for — U < f < 0, but the overall 
magnitude of the current decreases largely with the increasing 
interaction strength. 



purpose we first find (again assuming r Q (e) constant) 



n a (r) = 



dE- 



|A Q | y/E*-\A a \2 

K (\TA a \)-K (([3-\T\)\A a \), 



e -E\r\ _ e -E(p-\r\) 



(52) 



where Kq (x) is the modified Bessel function of the second 
kind. The expression for the critical current thus reads 



Ir, = — 



T L T R \A L A B \ 



r-0 r-0 

dry / dr 2 / dr^Hhiji ~ t 2 ) 
o Jo Jo 



X H R {T 3 )B(T 1 ,T2,T 3 )V{T 1 ,T 2 ,n). 



(53) 



We evaluated I c numerically for a number of different 
cases with the qualitatively same results showing that the 
coupling to oscillator mode(s) suppresses the magnitude 
of the Josephson current. There is no apparent differ- 
ence between the single mode and Ohmic heat bath case 
which is in a clear contrast with nonequilibrium trans- 
port studies where the character of the phonon spectrum 
plays a crucial role in the current- voltage characteristics. 
Below, just for simplicity, we only present results for the 
symmetric case |A^| = \A R \ = A at zero temperature 
T = and for infinite Coulomb interaction U — » oo. 
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FIG. 3: (Color online) Josephson current dependence on the 
temperature for U — > oo. We show the temperature depen- 
dence of the critical Josephson current for |£| < A for three 
different temperatures T = (/3 -> oo) (full line), /3A = 40 
(dash-dotted line), and /3A = 10 (dashed line with asterisks). 
In the last case which illustrates roughly the highest temper- 
ature achievable within the approximations used the analytic 
and numerical results are compared. 
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FIG. 4: (Color online) Josephson current for low-frequency 
phonons at T = 0, U — ► oo. As long as the spectrum of 
the phonon mode(s) is well below the superconducting gap, 
i.e., uo,c <C A, the critical Josephson current only depends 
on the integral £ c i = / °° du) J {oj) /uj = gu>o, c via the shifted 
dissipationless function iV(£ + e c \ sgn£) (dash-dotted line). 
Numerical results are shown for a single phonon mode with 
u>o ~ 0.5 A, g = 2 (dots), ujq — 0.1 A, g = 10 (asterisks), 
ujo = 0.05 A, g — 20 (pluses), and for the Ohmic bath with 
uj c — A, g = 1 (crosses). All of these cases yield basically the 
same result very well captured by the analytic expression. 



A. Low-frequency phonons 

If the spectrum of the oscillator mode(s) is well below 
the superconducting gap wo jC "C A we can find an ap- 
proximate analytic expression for the I C (S.) with the help 
of the function n(£) of the U — > oo case with no phonons, 
Eq. H48fl . To this end we study the current formula 153|) 
when we plug into it the expressions for B(t±, T2, t%) 141|) 
and 2?(ti,T2,T3) (|32|l and consider the above limit. For 
(<0we notice that the step functions of Eq. (|41|l and 
the fast decaying functions 7i Q (r) in Eq. (|53|) limit the 
relevant contributions to the three-dimensional integral 
to values of n, t 2 , r 3 small compared to min{l/|£|, 1/A}. 
For that reason one can perform the Taylor expansion in 
Eq. of the h(r) functions given by Eq. (|31J which 
gives 

%)«-|r| r dco ^ = -e cl \r\, (54) 
Jo u 

where we have defined the quantity 

u V V 

being the classical displacement energy when the oscil- 
lators are displaced by l v due to the force generated 
by a single excess electron. Putting the expansion into 
Eqs. (321 and (jHSJ we get 

r L r^A 2 e-n f r> p 

X H L (n - T 2 )H R (T 3 ) e (e-ci)(Ti-r 3 +r 2 ) ) 

(56) 

valid for for £ < 0. This in fact corresponds to the dis- 
sipationless case with the replacement £ — > £ — e c \. In 
particular, in the zero temperature limit we get I c (0 = 
-(rLr fl /7T 2 A)n[(£-e cl )/A] fore; < 0. For £ > we have 
to combine the vanishing prefactor exp(— /3£) with the di- 
vergent integral (for (3 — » oo) to get a finite result. This 
can be done by the substitutions r t - = /3 — Tj, i = 1, 2, 3 
which make the integrand relevant only for small val- 
ues of t/'s and, analogously to the previous derivation, 
one finally finds in the zero temperature limit, that 
7 C (0 = (2r L r fl /7r 2 A)n[(£ + e cl )/A] valid for £ > 0. In 
total, the critical Josephson current in case of coupling 
to low-frequency oscillator modes can be expressed as 
(remember the relation (|4(j[) between N(tf) and n(x)) 

^) = ^A^ + e cl sgn£). (57) 

This result is illustrated in detail in Fig. 0] 

We point out that even in case of finite temperature 
there is no residual temperature dependence due to cou- 
pling to the bath even though the low-frequency modes 
could be significantly populated. This is seen from the ex- 
pansion (|54|l where the thermal occupation factors ub{oj) 
cancelled out. Thus, the only temperature dependence 
would be the one stemming from the occupation factors 
exactly as in the dissipationless case, Eq. Q49f) and Fig. [3] 
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FIG. 5: (Color online) Josephson current for high-frequency 
single oscillator at T = 0, U — > oo,ljo = 20 A. For large 
enough ojo S> A, £ the Josephson current is expected to be 
suppressed just by the factor exp(— 2g) compared to the case 
without the phonon. This is shown in the figure where the 
numerical results for several values of the coupling constant 
g — 0.1 (asterisks), g = 0.5 (pluses), g = 1 (crosses), and 
g — 2 (dots) are compared with the analytic expressions given 
by the corresponding lines. 



B. High-frequency single oscillator mode 

For the high-frequency single phonon mode uiq A, £ 
we expect suppression of the critical Josephson current 
due to the fact that only transport through the ground 
oscillator state is allowed since the (virtual) involvement 
of the excited states would be further suppressed by a fac- 
tor A/nuiQ <C 1 (n = 1, 2, 3, . . .). The transport through 
the ground state is then diminished by the overlap fac- 
tors /oo = (0|e _jioPo |0)o, more precisely the Josephson 
current should be suppressed by a factor of order 




At T — and for ljq ^> A, £ the function h osc (r) — 
(?(e™" ' T ' — 1) changes very fast for small r's of the order 
r ~ 1/ujq which is irrelevant for the integral (|53l) . For 
larger r's the function h osc (r) ~ — g is constant and thus 
we get from Eq. I|32|) that the total effect of the phonon 
on the Josephson current is just a constant factor of e~ 2s 
multiplying the dissipationless case. This effect is shown 
in Fig.|3for different values of the coupling constant g at 
fixed wo = 20 A. The results do not depend on the value 
of luq provided it is high enough, which depends on the 
value of g, as expected (not shown). One can see that the 
approximation gets worse with increased g (for fixed u>q) 
which can be explained by the increased contribution of 
higher order virtual processes favored by the larger value 
of the coupling constant. 



FIG. 6: (Color online) Josephson current for high-frequency- 
cutoff Ohmic bath at T — 0, U — » oo and for several values 
of the coupling constant g and cutoff frequency w c . We see 
a significant suppression of the critical current already for 
intermediate coupling g — 1. 

C. Ohmic bath with high-frequency cutoff 

The case of Ohmic bath with large cutoff energy may 
be the physically most relevant one. Unfortunately, there 
is no simple semianalytic theory for this case and thus 
we have to rely mainly on the numerical results which 
are summarized in Fig. [BJ One should notice the very 
strong dependence on the coupling constant g. Even for 
intermediate coupling strength g ~ 1 we get a significant 
suppression of I c . 

We can give a qualitative explanation of the strong 
suppression due to the large frequency part of the phonon 
spectrum, uu > A. Roughly, we capture the effect of these 
high-frequency modes by a suppression factor similar to 
Eq. H58|) . This factor is estimated as 




Thus for large lo c and/or g the supercurrent is suppressed 
quite severely, in qualitative agreement with the results 
presented in Fig. |H1 

VI. SUMMARY AND DISCUSSION 

We have calculated the Josephson current through a 
single correlated level with coupling to external bosonic 
degrees of freedom representing, e.g., a system consist- 
ing of a molecular transistor with a number of internal 
vibrational modes and coupled to the phonons of the sub- 
strate. 

First, we have studied the case without coupling to vi- 
brations. This situation has been studied previously in 
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a number of papers, but we have derived new analytic 
formulae. The effect of the Coulomb interaction is to 
strongly suppress the Josephson current at the charge de- 
generacy point, and since this is where the junction also 
crosses over to 7r junction behavior, the critical current 
in fact goes to zero at this point. One could check this 
behavior using, e.g., nanotube devices coupled to super- 
conductors as in Refs. [ijj [l|J by tuning the gate voltage 
across the charge degeneracy point. Also the tempera- 
ture dependence predicted here could be experimentally 
verified. 

In the second part of the paper we have included the 
coupling to environmental modes and discussed different 
limits. Coupling to low-frequency phonons does not have 
a severe influence on the Josephson coupling, since it only 
shifts the argument of the dissipationless formula from 
the single particle energy by the classical displacement 
energy. In contrast, a strong coupling to high-energy 
phonons suppresses the supercurrent quite substantially. 
This is because the oscillators are displaced twice during 
the transfer of Cooper pair, and because the supercurrent 
requires the final state of the oscillator to be identical to 
the initial one, the transfer is suppressed by twice the ex- 
ponential Franck-Condon factor e~ 2g . Therefore it might 
be difficult to observe Josephson tunneling current in de- 
vices with a strong electron- vibron coupling. 

Throughout the paper we have used lowest-order per- 
turbation theory in the tunneling coupling. For stronger 
tunneling coupling one expects the correlation effect due 
to the vibrations to become smaller, because the charge 
on the level is no longer well defined. This means that 
a mean-field treatment becomes adequate, when T > e c \. 
How to describe this transition theoretically is an inter- 
esting problem. 
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APPENDIX A: NONINTERACTING CASE U = 

The noninteracting case is exactly solvable for any val- 
ues of parametersiS&Si Here, we only give a brief sketch of 
the solution and the summary of the results in the limit 
of small r relevant for our study. Since the model for 
U = has a quadratic Hamiltonian, it can be solved 
by the equation of motion technique for the Matsub- 
ara Green function. We introduce an infinite vector 

a = (dT'4' c feiT' c -fci|' c fc«T' c -fei?i) T generalizing the 
standard Nambu formalism to the case of the dot plus 
two leads. Defining the corresponding thermal Green 
function Q(t) — —(T T (a(r)cv)) satisfying the equation 



of motion 4-Q{r) = -S(r)l + M- Q(t) with the matrix 



iM = 



/ 







— t* 
l kL 





l kR 


\ 










t-kL 





t-kR 




-tkL 





~£kL 


Al 













t* 

b -kL 


Al 


£kL 










-tkR 











~£kR 


Ar 


V 





t* 

b -kR 








A* 

^R 


£kR ) 



(Al) 



we can express the Josephson current as 
II = -4Im^i_ fcL £/4 2 (T -> 



0+) 



(A2) 



(factor of 2 for spin degeneracy). Going to the frequency 
picture (-4- — > — iui n ) and using the partitioning scheme 



A c 
d B 



{A-cB^d)- 1 
-B- l d{A-cB- 1 ^- 1 



-A _1 c(fl - dA^c)- 1 
d)- 1 (B-dA^c)- 1 

together with the wide-band approximation r Q = const, 
and assuming the symmetric case T l = Tr = T, = 
Ae 1 *, Ar = A this set of linear equations gives for the 
Josephson current 



/,=Imi£ e — 



A 2 )D(ui n ) 



(A4) 



with 



D{w n )=ul{l + 



A 2 r 2 cos 2 (0/2) 
(c 2 +A 2 ) ■ 



(AS) 

The zeros of -D(cOn) determine the Andreev bound states 
discussed in Ref . |3- In the limit T — > which we consider 
here the lowest order contribution to II is proportional 
to r 2 and can be obtained by setting T = in D(uj n ). 
The sum can then be easily performed yielding 



77, = — sin ( 



/tanh(/?£/2) tanh(/3A/2) 



2(l-(£/A) 2 )V € 



(A6) 

For UbT <C A this expression is proportional to V 2 /UbT 
at £ = 0. This would diverge in the T = limit. However, 
this divergence is just an artifact of our perturbation the- 
ory and the exact evaluation of the full expression (|A4|) 
would give a finite result even at T = 0, £ — > 0. The exact 
result is essentially identical to our approximate one un- 
less £, ksT < r when the exact result gives a saturation 
of the maximum Josephson currenfcS&Sil J c s» T around 
£ = for r <C A. For nonsymmctric coupling to the 
leads, i.e., T L ^ T R but still |A L | = |A H | = A, we get 
within the discussed lowest order approximation in T's 
the same result just with the replacement T 2 — > r^r^ 
which is used in the main text, Sec. lTVAl 



APPENDIX B: FINITE U, EVALUATION OF n(x) 

In this appendix we calculate the critical Josephson 
current for finite value of the Coulomb interaction U 
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in the limit of zero temperature T — and symmetric 
gap | A^l = |A#| = A. The critical current is given by 
Eq. 142(1 with the function N{£) being determined in this 



case by Eqs. (|43fl . (|50|) . and (|51|l which further simplify in 
the considered limit of zero temperature and symmetric 
gap so that we can write for N(£): 



N(0 = 



f OO pOO 

dE / dE'- 



2A 2 



A VE 2 - AVE' 2 - A 2 
arccos 2 (£/A) 



1 



OO POO 

dE / dE' 



(£ + [//2)(l-(£/A) 2 ) 
A 2 



(E + E') {E + (E' + (2£ + + + #') 
for £>0, 

1 1 



(Bla) 



(E + E')V# 2 - AVE' 2 - A 2 



{E-£)(E'-£) {E + t + U){E' + i + U) 



{E-i){E + i + U) 



OO /" OO 

dE / dE' 



i(-£/A) + n((£ + CO/A) + 4n 2 (-£/A, + f/)/A) 



for - [7 < £ < 0, 



(Bib) 



A 



2A 2 



VE 2 - AVE' 2 - A 2 



(E + E') + E) (£ + [/ - E>) (2£ + U){£ + U- E)(£ + U-E') 



= 2 



'n(|£ + ET|/A) 



arccos 2 (|£ + [/|/A) 
(e + C//2)(l-(|e + C/|/A) 2 ) 



for £ < -C7. 



(Blc) 



Here, n(ir) is defined by Eq. I|47(l and 

, , r°° du r°° dv 



i y/u 2 - 1 k Vv 2 - 1 (u + v)(u + x){u + 
hi(u + Vu 2 - 1) 



fit/ 



(u 2 — 1)(« + x)(u + y) 



(B2) 



Again, the appropriate analytic continuation of arccos x 
for x > 1 is understood (see discussion below Eq. iJUJl). 
While neither of the integral definitions of n(x) (Eq. (|47(l ^ 



or n 2 (x, y) (Eq. l|B2|ll seems to lead to an explicit formula 
we still could, however, use the above equations for in- 

direct evaluation of n(x). In particular, since the first 

y) equation (|Bla|) is valid for £ > for any value of U 
including both the limits C7 — * and U — > oo we can 
use the former limit for the evaluation of the latter one. 
The noninteracting case U = is described by Eq. 140|) 
and, thus, we come to the simple result (|48|l (in the limit 
f3 — > oo). The function 712(2?, 2/) can be easily calculated 
by numerical evaluation of the integral in Eq. 1|B2(I for 
x,y > -1. 
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